#include <iostream>
#include <fstream>
#include "optical.h"
#include "density_generator.h"
#include "full_sed_grid.h"
#include "random.h"
#include "boost/program_options.hpp"
#include "scatterer.h"
#include "dust_model.h"
#include "full_sed_emergence.h"
#include "montecarlo.h"
#include "vecops.h"
#include "xfer.h"
#include "emission-fits.h"
#include "emission.h"
#include "mcrx.h"
#include "wd01_Brent_PAH_grain_model.h"
#include "shoot.h"
#include "imops.h"

#include "mpi.h"
#include "arepo_grid.h"

using namespace std;
using namespace mcrx;
using namespace boost;
namespace po = boost::program_options;


typedef local_random T_rng_policy;
typedef scatterer<polychromatic_scatterer_policy, T_rng_policy> T_scatterer;
typedef T_scatterer::T_lambda T_lambda;
typedef T_scatterer::T_biaser T_biaser;
typedef dust_model<T_scatterer, cumulative_sampling, T_rng_policy> T_dust_model;
typedef full_sed_emergence T_emergence;
typedef emission<polychromatic_policy, T_rng_policy> T_emission;
typedef full_sed_grid<arepo_grid> T_grid;
typedef T_grid::T_cell_tracker T_tracker;


int main(int argc, char** argv)
{
  // Declare the supported options.
  po::options_description desc("Optically thick cloud with central point source test problem. Allowed options");
  desc.add_options()
    ("output-file", po::value<string>(), "output file name")
    ("n_rays", po::value<long>()->default_value(1000000), "number of rays")
    ("n_rays_temp", po::value<long>()->default_value(1000000), 
     "number of rays for dust temp iteration")
    ("n_rays_ir", po::value<long>()->default_value(1000000), 
     "number of rays for final dust emission")
    ("n_threads", po::value<int>()->default_value(12), "number of threads")
    ("maxlevel", po::value<int>()->default_value(6), 
     "max grid refinement level")
    ("n_ir_forced", po::value<int>()->default_value(0), 
     "number of forced scatterings for IR emission")
    ("vis_ref", po::value<T_float>()->default_value(0.9e-6),
     "stellar radiation reference wavelength")
    ("ir_ref", po::value<T_float>()->default_value(40e-6),
     "dust emission reference wavelength")
    ("i_rr", po::value<T_float>()->default_value(1e-2),
     "intensity where Russian Roulette starts")
    ("rel_tol", po::value<T_float>()->default_value(1e-3), 
     "relative ir equilibrium tolerance")
    ("abs_tol", po::value<T_float>()->default_value(1e-3), 
     "absolute ir equilibrium tolerance")
    ("source_temp", po::value<T_float>()->default_value(5800),
     "temperature of central blackbody source")
    ("source_lum", po::value<T_float>()->default_value(1),
     "luminosity of central source (in Lsun)")
    ("radius", po::value<T_float>()->default_value(5),
     "radius of cloud (in AU)")
    ("tau", po::value<T_float>()->default_value(100),
     "V-band optical depth to center of cloud")
    ("tau_tol", po::value<T_float>()->default_value(100),
     "cell optical depth refinement tolerance")
    ("seed", po::value<int>()->default_value(42), "random number seed")
    ("use_blackbody", po::value<bool>()->default_value(false), "Use ideal blackbody grains instead of a physical dust model ")
    ("grain_data_dir", po::value<string>()->default_value("/home/patrik/dust_data/crosssections"), "grain data directory")
    ("help", "print this page")
    ;

  po::positional_options_description p;
  p.add("output-file", 1);

  po::variables_map opt;
  po::store(po::command_line_parser(argc, argv).
	    options(desc).positional(p).run(), opt);
  po::notify(opt);    

  if (opt.count("help")) {
    cout << desc << "\n";
    return 1;
  }
  if(!opt.count("output-file")) {
    cout << "Must specify output file name on command line" << endl;
    return 1;
  }
  
  cout << "Running ir_thick_test with parameters:\n";
  cout << "output-file = " << opt["output-file"].as<string>() << '\n'
       << "n_rays = " << opt["n_rays"].as<long>() << '\n'
       << "n_rays_temp = " << opt["n_rays_temp"].as<long>() << '\n'
       << "n_rays_ir = " << opt["n_rays_ir"].as<long>() << '\n'
       << "n_threads = " << opt["n_threads"].as<int>() << '\n'
       << "maxlevel = " << opt["maxlevel"].as<int>() << '\n'
       << "n_ir_forced = " << opt["n_ir_forced"].as<int>() << '\n'
       << "vis_ref = " << opt["vis_ref"].as<T_float>() << '\n'
       << "ir_ref = " << opt["ir_ref"].as<T_float>() << '\n'
       << "i_rr = " << opt["i_rr"].as<T_float>() << '\n'
       << "rel_tol = " << opt["rel_tol"].as<T_float>() << '\n'
       << "abs_tol = " << opt["abs_tol"].as<T_float>() << '\n'
       << "source_temp = " << opt["source_temp"].as<T_float>() << '\n'
       << "source_lum = " << opt["source_lum"].as<T_float>() << '\n'
       << "radius = " << opt["radius"].as<T_float>() << '\n'
       << "tau = " << opt["tau"].as<T_float>() << '\n'
       << "tau_tol = " << opt["tau_tol"].as<T_float>() << '\n'
       << "seed = " << opt["seed"].as<int>() << '\n'
       << "use_blackbody = " << opt["use_blackbody"].as<bool>() << "\n\n";

  const int n_threads = opt["n_threads"].as<int>();
  const long n_rays =opt["n_rays"].as<long>();
  const long n_rays_ir =opt["n_rays_ir"].as<long>();
  const long n_rays_temp = opt["n_rays_temp"].as<long>();
  const int grid_n = 10;
  const int maxlevel= opt["maxlevel"].as<int>();
  const int n_forced=1;
  const int n_forced_ir=opt["n_ir_forced"].as<int>();
  const T_float vis_ref_lambda = opt["vis_ref"].as<T_float>();
  const T_float ir_ref_lambda = opt["ir_ref"].as<T_float>();
  const T_float rel_tol =opt["rel_tol"].as<T_float>();
  const T_float abs_tol =opt["abs_tol"].as<T_float>();
  const T_float tau = opt["tau"].as<T_float>();
  const T_float tau_tol = opt["tau_tol"].as<T_float>();
  const string output_file(opt["output-file"].as<string>());
  const T_float i_rr=1e-2;
  const T_float cameradist=1e4;
  const T_float source_temp = opt["source_temp"].as<T_float>();
  const T_float source_lum = opt["source_lum"].as<T_float>();

  Preferences prefs;
  prefs.setValue("mcrx_data_dir", string("~/sunrise/src/"));
  prefs.setValue("grain_data_directory", string("~/data/dust_data/crosssections"));
  prefs.setValue("wd01_parameter_set", string("DL07_MW3.1_60"));
  prefs.setValue("template_pah_fraction", 0.5);
  //  prefs.setValue("wd01_parameter_set", string("LMCavg_20"));3

  prefs.setValue("use_dl07_opacities", true);
  //prefs.setValue("use_dl07_opacities", false);
  prefs.setValue("n_threads", n_threads);
  prefs.setValue("use_grain_temp_lookup", true);
  prefs.setValue("use_cuda", false);

  if (prefs.defined("mcrx_data_dir"))
    mcrx::particle_base::inp_dir = 
      word_expand (prefs.getValue("mcrx_data_dir", string ())) [0];

  mcrx::seed(opt["seed"].as<int>());
  T_unit_map units;
  units ["length"] = "kpc";
  units ["mass"] = "Msun";
  units ["luminosity"] = "W";
  units ["wavelength"] = "m";
  units ["L_lambda"] = "W/m";

  CCfits::FITS f("grid_010_arepo.fits", Read);

  array_1 lambda;
  ExtHDU& lambda_hdu = open_HDU (f, "LAMBDA");
  read(lambda_hdu.column("lambda"), lambda);

  // load grain model
  T_dust_model::T_scatterer_vector sv;
  sv.push_back(boost::shared_ptr<T_scatterer> 
	       (new wd01_Brent_PAH_grain_model<polychromatic_scatterer_policy, mcrx_rng_policy>(prefs, units)));

  T_dust_model model (sv);
  model.set_wavelength(lambda);
  vector<T_float> lambdas(lambda.begin(), lambda.end());

  const int lambda550 = lower_bound(lambdas.begin(), lambdas.end(),
				    0.551e-6) - lambdas.begin()-1;
  // find opacity normalization for whatever dust we're using
  const T_float unit_opacity_rho = 1.0/sv[0]->opacity()(lambda550);

  // emission
  auto_ptr<T_emission> em;
  em.reset(new full_sed_particles_emission<cumulative_sampling, local_random>(open_HDU(f, "PARTICLEDATA"), false ));


  cout << "setting up emergence" << endl;
  vec3d translate_origin;
  ExtHDU& shdu=open_HDU(f,"SFRHIST");
  shdu.readKey("galcenter1X",translate_origin[0]);
  shdu.readKey("galcenter1Y",translate_origin[1]);
  shdu.readKey("galcenter1Z",translate_origin[2]);
  cout << "Translating camera origin to " << translate_origin << '\n';
  std::vector<std::pair<T_float, T_float> > cam_pos;
  cam_pos.push_back(std:: make_pair (0.1, 0.01) );
  cam_pos.push_back(std:: make_pair (60*constants::pi/180, 0.01) );
  cam_pos.push_back(std:: make_pair (77*constants::pi/180, 0.01) );
  cam_pos.push_back(std:: make_pair (90*constants::pi/180, 0.01) );
  cam_pos.push_back(std:: make_pair (90*constants::pi/180, 89*constants::pi/180) );
  auto_ptr<T_emergence> e(new T_emergence(cam_pos, cameradist, 30, 300, 
					  translate_origin));

  std::vector<T_rng::T_state> states= generate_random_states (n_threads);
  T_rng_policy rng;

  MPI_INIT(&argc, &argv);
  arepo::load_snapshot("snapshot_010","parfile-sunrise",n_threads);

  shared_ptr<arepo_grid<absorber<array_1> > > gg(new arepo_grid<absorber<array_1> >());
  uniform_density_generator dg(0.01,0.4);
  T_grid gr(gg, dg, 25);
  gr.reset();

  xfer<T_dust_model, polychromatic_grid<arepo_grid>, T_rng_policy> x (gr, *em, *e, model, rng, i_rr, 0,n_forced,10.);
  long n_rays_actual = 0;
  dummy_terminator t;
  const int vis_ref = lower_bound(lambdas.begin(), lambdas.end(),
				  vis_ref_lambda) - lambdas.begin()-1;
  T_biaser b(vis_ref);
  shoot (x, scatter_shooter<T_biaser> (b), t, states, n_rays,
  	 n_rays_actual, n_threads,0, "");

  T_float normalization = 1./n_rays_actual;

  {
    FITS output("!"+opt["output-file"].as<string>(), Write);
    e->write_images(output,normalization,units, "-STAR", false, false);

    // calculate integrated SEDs
    Table* iq_hdu=output.addTable("INTEGRATED_QUANTITIES",0);
    iq_hdu->addColumn(Tdouble, "lambda", 1,
		      units.get("L_lambda"));
    write(iq_hdu->column("lambda"),lambda,1);
    for (int i=0; i<cam_pos.size(); ++i) {
      integrate_image(output,
		      "CAMERA"+boost::lexical_cast<string>(i)+"-STAR",
		      4*constants::pi*cameradist*cameradist,
		      "L_lambda"+boost::lexical_cast<string>(i),
		      lambda,
		      "",
		      units,
		      dummy_terminator());
    }
  }

}
